Absolute excited state molecular geometries revealed by resonance Raman signals

Ultrafast reactions activated by light absorption are governed by multidimensional excited-state (ES) potential energy surfaces (PESs), which describe how the molecular potential varies with the nuclear coordinates. ES PESs ad-hoc displaced with respect to the ground state can drive subtle structural rearrangements, accompanying molecular biological activity and regulating physical/chemical properties. Such displacements are encoded in the Franck-Condon overlap integrals, which in turn determine the resonant Raman response. Conventional spectroscopic approaches only access their absolute value, and hence cannot determine the sense of ES displacements. Here, we introduce a two-color broadband impulsive Raman experimental scheme, to directly measure complex Raman excitation profiles along desired normal modes. The key to achieve this task is in the signal linear dependence on the Frank-Condon overlaps, brought about by non-degenerate resonant probe and off-resonant pump pulses, which ultimately enables time-domain sensitivity to the phase of the stimulated vibrational coherences. Our results provide the tool to determine the magnitude and the sensed direction of ES displacements, unambiguously relating them to the ground state eigenvectors reference frame.

Photochemical and photophysical reactions are ruled by excited state (ES) potential energy surfaces (PESs), which drive nuclear motions and initiate the molecular rearrangements along specific vibrational degrees of freedom. The nuclei respond to photo-excitation altering their equilibrium positions, i.e., the molecule undergoes a geometrical rearrangement which, in turn, modifies the physical and chemical properties of the system. Nature has tailored electronically excited PESs in which the vibrational structures are specifically modified from the ground state (GS) equilibrium configuration to efficiently convert the absorbed light energy into molecular rearrangements, driving the biological function and the photo-chemistry by bond length modifications, torsional re-orientations, formation or rupture of chemical bonds. These can be rationalized by the multidimensional displacement D between PESs, i.e., the vector identifying the new equilibrium position in the ES projected onto the eigenvectors (Q j ) of the GS with respect to its minimum. Large absolute values of specific displacement components (d j = D ⋅ Q j , commonly referred to as displacements) identify the vibrational eigenmodes which are more involved in the excitation and eventually their couplings to the specific electronically excited state. Most importantly, the sensed direction of D determines whether the changes in the equilibrium bond distances and angles move closer or away two functional groups, ruling the physicalchemical properties on the ES and, thus, its knowledge is of utmost importance.
For example, while the ultrafast olefinic photoisomerization of stilbene is determined by an increase of the central C=C bond length with out-of-plane motion of the two ethylenic hydrogens 1 , the mechanism of vision is habilitated by a sub-picosecond retinal isomerization of the rhodopsin molecule, initially driven by positive and negative changes of different dihedral angles 2 . In diaryl thiophenes, a different ES dynamics is underpinned by an increase of the distance between the carbon (connected to phenyl) and the sulfur atoms in the thiophene as opposed to a decrease of the dihedral angle between phenyl and thiophene rings 3,4 . The displacement sense also plays a role in the excited states of conjugated push-pull systems, where the alternation of expansion and contraction in specific bond lengths is suggestive of an enhanced zwitterionic character 5 . Finally, the design of molecular motors requires an asymmetry in the excited state PES along the dihedral angle with respect to the ground state, which can be provided, for instance, by chirality, in order to achieve a preferential direction and obtain photoinduced molecular rotary motion 6,7 . Measuring excited state displacements is also important to validate quantum chemistry models for nonlinear optical effects in molecular compounds. These examples demonstrate the importance of determining the sensed direction of the displacements within the Franck Condon (FC) description of the excited state. This information, however, is highly convoluted and difficult to access from the static spectra obtained by absorption or other linear techniques. In particular, reorientation along opposite senses often leads to very similar experimental signatures, motivating a quest for suitable methods to access this quantity.
Spontaneous Raman spectroscopy represents a powerful technique to study the vibrational spectrum of molecular and solid state systems in the GS. Tuning the excitation wavelength to match the energy of an allowed electronic transition to a given ES enables resonant Raman (RR) spectroscopy [8][9][10][11][12][13][14][15][16] and results in a cross section proportional to the square modulus of the product between the electronic transition moments and the vibrational GS-ES overlaps 11,[17][18][19][20][21][22] , ultimately related to their relative PES geometries. This fundamentally limits the RR sensitivity uniquely to the absolute values 23 of the ES displacements 24,25 . The missing sign information hampers elucidating the molecular reorientations at the basis of photo-physical and photochemical phenomena 10 . In addition, measuring and isolating the RR response under identical conditions for a sequence of narrowband experiments is often a challenging task, further hampered by the fluorescence emission, which typically overwhelms the weak spontaneous Raman signatures.
Nonlinear time-domain Raman schemes, in particular impulsive stimulated Raman scattering (ISRS) 1,26-31 , offer the chance to circumvent some of the spontaneous Raman limitations. ISRS employs a pair of temporally well separated pulses to investigate vibrational matter properties: two interactions with a femtosecond Raman pump (RP) generate vibrational coherences, triggering nuclear wavepacket motions and modulating the sample transmissivity, which is then measured by a broadband probe pulse (PP) as a function of its temporal delay ΔT with respect to the RP. This technique records the Raman response directly in the time-domain, thereby simultaneously accessing the amplitude and the phase of the vibrational oscillations [32][33][34][35][36][37][38][39][40][41][42][43] . Fast Fourier transforming (FFT) the signal with respect to ΔT recovers the conventional frequency-domain Raman spectrum. In view of the heterodyne detection, the measured signal is generated on top of the PP and hence the incoherent isotropic fluorescent background is efficiently suppressed 20,44 . Moreover, the time-domain acquisition, performed with temporally separated pump and probe pulses, ensures a Raman response free from nonlinear artifacts caused by cross-phase modulation 30,45,46 , which typically affects frequencydomain stimulated Raman approaches 46-48 . In this work, we introduce an experimental protocol, based on a two-color broadband impulsive stimulated Raman scattering (ISRS) scheme, that can determine the real and imaginary parts of the nonlinear RR susceptibility, revealing absolute ES geometries.

Results and discussion
Design of a two-color ISRS scheme to measure the sense of ES displacements To take advantage of the phase sensitivity of the ISRS response, we analyze the signal in the time-domain before Fourier transforming, since evaluating the squared absolute value of the FFT spectra loses the phase information contained in the time-domain data. We employ an off-resonant pump to induce vibrational coherences only on the initially populated electronic level and a probe pulse resonant with a targeted electronic transition. Importantly, the off-resonant pump is critical in that it ensures a stimulation of the vibrational coherences via the molecular polarizability, and not via the Franck-Condon overlap integrals, which are probed solely by the PP. Hence, under the adopted scheme, the ISRS signals are proportional, through the off-resonant GS polarizability, to the amplitude of the GS-ES Franck-Condon overlap integrals and not to their square modulus 49 , as opposed to spontaneous Raman spectroscopy. This approach makes the RR ISRS response, monitored as a function of the PP wavelength, sensitive to the sign of PES displacements. The concept is demonstrated in Fig. 1 for the illustrative case of a diatomic molecule, where the unique degree of freedom is the distance between the nuclei, which can be either increased or decreased on the ES (panel a). Importantly, ES molecular geometries corresponding to an increase or a decrease of the nuclear bond-length have FC overlaps with the same magnitude and opposite sign. Hence, the ISRS signal, recorded with the two-color non degenerate scheme presented in this work (panel b), is antisymmetric with respect to the GSto-ES atomic distance modification (panels c-d) and directly reveals the ES displacement. It is worth to stress that, in principle, the assignment of the sign for a given ground state eigenvector Q, obtained by diagonalizing the Hessian matrix, is arbitrary in the sense that each arrow identifying the direction of an atomic vibration may be reversed (Q 0 = À Q, panel e). Importantly, since the proposed ISRS scheme accesses the product between the GS polarizability derivative ∂α ∂Q and the FC factors, reversing the sign for a given ground state eigenvector changes both the sign of the polarizability derivative as well as the sign of the FC overlaps (panels f-g), making the experimental signals invariant on the selection of such reference frame. Hence, at variance with the spontaneous Raman approach, the sign of the ES displacements relative to the GS eigenvectors can be defined and the ES equilibrium configuration along a given vibrational mode is uniquely identified from the measured ISRS Raman excitation profiles (REPs).
The broadband ISRS response can be evaluated -in a semiclassical framework-through a perturbative expansion of the molecular density matrix in powers of the classical electric fields EðtÞ = P i E i ðtÞe Àiω i t + c:c:, using a diagrammatic approach to visualize the concurring processes that generate the time-domain Raman signal 44,50,51 . Particularly, since the RP can act both on the ket and on the bra side of the density matrix, two classes (A/B) of Feynman diagrams (reported in Fig. 2a, b) have to be considered to calculate the third-order polarizabilities, P ð3Þ A ðω, ΔTÞ and P ð3Þ B ðω, ΔTÞ respectively, that originate the system response 51,52 . For small modifications of the PP across the sample, the spectrally dispersed ISRS signal S(ω, ΔT), defined as the normalized difference between the transmitted PP spectrum in presence (I RP On P ðω, ΔTÞ) and in absence (I P (ω)) of the RP, can be expressed as the imaginary part (ℑ) of the ratio between the total third-order polarization P ð3Þ ðω, ΔTÞ = P ð3Þ A ðω, ΔTÞ + P ð3Þ B ðω, ΔTÞ and the probe field spectral envelope E P (ω), accordingly to Sðω, ΔTÞ = I RP On P ðω, ΔTÞ À I P ðωÞ I P ðωÞ / À= P ð3Þ ðω, ΔTÞ=E P ðωÞ It is worth to stress that, as the RR process stimulated by the probe pulse projects the vibrational wavefunction onto the ES manifold, both the P ð3Þ A ðω, ΔTÞ and the P ð3Þ B ðω, ΔTÞ terms contain sumover-states expressions 44,53 , and are therefore sensitive to the ES geometry.
For an electronically off-resonant RP, the radiation-matter interaction Hamiltonian involved in the preparation of the vibrational coherences is 54 H ðRPÞ I = Àα Á jE R j 2 , where α indicates the molecular polarizability, while for the resonant PP detection it reads as H ðPPÞ I = À μ Á E P , where μ is the dipole operator. Indicating the frequency difference between the j-i levels as ω ji = ω j − ω i , the third order polarizations related to the P ð3Þ A ðω, ΔTÞ and P ð3Þ B ðω, ΔTÞ terms can be written as where is the vibrational coherence preparation function (that is the same for the two A/B processes for an off-resonant pump, as detailed in Supplementary Note 1), which depends on the Raman field E 0 R ðωÞ and on the molecular polarizability via ∂α ∂Q g 0 . The absolute value of the latter can be conveniently obtained by an off-resonant Raman experiment, while its sign encodes the structure of the GS eigenvectors which, in turn, define the normal modes reference frame. Importantly, the information regarding the excited state displacement (along the considered ω g 0 g normal mode) is fully contained in the R S g 0 ðωÞ and R AS g 0 ðωÞ REPs, which are probed respectively by the P ð3Þ A and the P ð3Þ B pathways. Since the R AS g 0 term is shifted by one vibrational quantum with respect to R S g 0 , they corresponds to the anti-Stokes and Stokes complex Raman excitation profiles: where the summation over k in Eq. (4) takes into account the ES vibrational manifold andω ji = ω ji À iΓ ji , which depends on the dephasing rate Γ ij = ðT ji vib Þ À1 of the j i h j coherence. It is worth to stress that, while the A/B diagrams reported in Fig. 2a, b do not correspond to the spontaneous Raman Stokes/anti-Stokes processes, the square moduli of the REP reported in Eq. (4) represent the information that can be accessed by spontaneous Raman spectroscopy 11 . Crucially, R S g 0 ðωÞ and R AS g 0 ðωÞ include the product of two dipole matrix elements, namely μ g 0 e k μ e k g . As shown in Fig. 2d, this product is an odd function of the displacement between the ground and the excited PESs along the considered normal mode: hence, accessing the sign of the real/ imaginary part of the REPs, which depend linearly on the μ g 0 e k μ e k g product (cf. Fig. 2e), provides the chance for unambiguously determining the sign of ES displacements along the considered normal modes.
Based on the odd symmetry of the complex valued Raman excitation profiles with respect to the displacement (Fig. 2e), the sign of this latter can be in fact identified by the ISRS signal. Aiming at directly accessing the real and imaginary REPs components, which in the measured signal are entangled in a linear combination (see Eq. (6) in the Methods section), an appropriate tuning of the probe chirp (C 2 , related to the group delay dispersion GDD as C 2 = GDD=2) can be directly exploited. Indeed, since the P ð3Þ A and P ð3Þ B terms originate from Since the ISRS response is proportional to the product between the ground state (GS) polarizability derivative and the Franck-Condon (FC) overlaps, the measured signals S(λ) uniquely determines the ES geometry and the sense of the geometrical rearrangement c, d, being antisymmetric with respect to the nuclear distance modification. Notably, while the sense of the GS eigenvector is arbitrary, as both outgoing (Q) or in-going arrows (Q 0 = À Q) can be selected e, reversing such a sense changes both the sign of the polarizability derivative as well as the sign of the FC overlaps f, g, making the experimental signals invariant on the selection of the reference frame. For example, if the complex REP reconstructed from the ISRS signal is the one reported in c, the consequences of the two possible eigenmode choices e are depicted in the corresponding f: the eigenmode direction Q is associated to a positive ∂α ∂Q and a positive displacement d, while the direction Q 0 = À Q corresponds to a negative polarizability derivative and a negative d. Crucially, the two formal descriptions yield to the same physical observable, i.e., an increase of the bond length in the excited state (highlighted in red). On the other hand, the measured REP in c excludes an ES bond-length shortening (blue), as the corresponding descriptions shown in g are both associated to the signal in d, which has a reversed sign with respect to the measured one.
interactions with different PP spectral components 52,55,56 , C 2 introduces a desired complex phase factor between the A and B pathways, providing the chance to directly measure the displacement sense and the real/imaginary parts of the Stokes/anti-Stokes complex REPs. As summarized in Table 1 and detailed in the Methods section, two values of the probe chirp, namely C 2 = m 2π T 2 vib and C 2 = 1 2π m + 1 4 À Á T 2 vib (with m = 0, ± 1, ± 2, . . . ), and two values of the dechirped temporal delaỹ T = n T vib andT = ðn + 1 4 Þ T vib (corresponding to ϕ = 2nπ and ϕ = 2ðn + 1 4 Þπ, respectively, with n = 0, 1, 2, . . . ) are required to obtain the real and the imaginary parts of R S g 0 ðωÞ and R AS g 0 ðωÞ.

Impulsive Raman experiments and reconstruction of complex REPs
As a road test for the capabilities of the presented approach, we applied the two-color broadband ISRS scheme to study the excited state PES of Rhodamine B (RhB), a staining fluorescent dye extensively used for chemical and biological applications. In view of its high fluorescence quantum yield 57 and small Stokes-shift 58 , the fluorescent background overwhelms the spontaneous Raman response, which hence can be efficiently measured only for pump excitation wavelengths far from the absorption maximum 59 .
The ISRS spectra of RhB dissolved in methanol have been measured at ambient temperature as a function of the PP wavelength λ P , for ΔT 2 À0:3, 3 ½ ps, with a 10 fs time-step; in Fig. 3a, b we report the timedomain traces recorded using a PP with a vanishing chirp C 2 ≈ 0. The RP is tuned to be electronically off-resonant (≈700 nm) and stimulates the normal modes with a not vanishing polarizability derivative, while the PP spectral profile covers the entire sample absorption (cf. Fig. 3e). By fast Fourier transforming over the RP-PP time delay, the frequencydomain spectrum (reported in Fig. 3c, d) is extracted, identifying a strong mode at 620 cm −1 and other weaker Raman bands at 490, 735, 1280 and 1360 cm −1 . Indeed, the~18 fs temporal duration of the pump pulse guarantees efficient stimulation of coherences up to 1400 cm −1 .
In order to reconstruct the time-domain vibrational response of the single Raman modes, the ISRS traces are further fine scanned around ΔT 2 0:35, 0:85 ½ ps with a 4 fs time-step, and then fitted as the sum of sinusoidal terms, by the equation Sðλ, ΔTÞ = P k A k ðλÞ sin ω g 0 k g ΔT + ϕ k ðλÞ h i . The frequencies ω g 0 k g are fixed parameters, determined from the peak positions of the FFT spectra, while the extracted amplitudes A k (λ) and phases ϕ k (λ) are reported in Fig. 4 for different probe chirps (C 2 ) as a function of λ = λ P . We note that for C 2 ≈ 0, all the modes exhibit a π phase shift around the absorption maximum at 550 nm, where the oscillation amplitude vanishes.
The fitted A k (λ) and ϕ k (λ) parameters were then used to reconstruct the time-domain ISRS response of the different vibrational modes. In Fig. 5a, b, this is evaluated for the 735 cm −1 mode, testifying the sensitivity of the scheme also to isolate weak Raman signals (see Fig. 3c, d). Here two PP chirp values are considered, namely C 2 = -4 fs 2 and C 2 = 71 fs 2 ≈ 1 8π T 2 vib , which closely correspond to the values reported in Table 1 (with m = 0) for separately accessing the real and imaginary parts of the mode complex REP. In panels c, d, the signal is then evaluated as a function of the dechirped temporal delayT, while four slices of such maps atT=T vib = n andT=T vib = ðn À 1=4Þ are reported in Table 1 | By properly tuning the RP-PP temporal delay ΔT and the probe chirp C 2 , it is possible to access separately the real and the imaginary part of the sum and the difference of the complex REPsT = n T vibT = ðn + 1 4 ÞT vib C 2 = m 2π T 2 vib <½R S g 0 ðωÞ À R AS g 0 ðωÞ =½R S g 0 ðωÞ + R AS g 0 ðωÞ , prepare a vibrational coherence on the ket a or on the bra side b of the density matrix. Interactions with the ket/bra side of the density matrix and the (classical) electromagnetic fields are presented by continuous/dashed vertical arrows, respectively. The sample density matrix is then promoted on the electronic excited state (ES) manifold e k via an interaction, proportional to FC overlaps μ ge k μ e k g 0 , with the broadband resonant probe pulse, enabling to record upon a free induction decay the ground state (GS) Raman response in the time-domain. Onedimensional projections of the ES PES along a specific normal mode are reported in c. Since the a/b pathways responsible for the measured signal depend on the product between two dipole moments μ ge k μ e k g 0 (reported in d for a single mode) and not on their absolute squared value, the experimental response, reported in arbitrary units in e, is sensitive to the linear combination of complex Stokes/anti-Stokes REPs (R S and R AS , respectively), which reveal the sign of the molecular displacement.
(e, f) for n = 12 as a function of the probe wavelength. Specifically, in the chosen reference frame (with ∂α ∂Q g 0 < 0, cfr. Supplementary Fig. 6), panel (e) shows À<½R S g 0 ðωÞ À R AS g 0 ðωÞ and =½R S g 0 ðωÞ + R AS g 0 ðωÞ as continuous green and yellow lines, respectively, while panel (f) =½R S g 0 ðωÞ À R AS g 0 ðωÞ and <½R S g 0 ðωÞ + R AS g 0 ðωÞ as blue and purple traces. To verify the reconstructed profiles accuracy, the experimental results are compared with the REPs modeled directly from the sample absorption spectrum by using the transform theory 60 (details are reported in the Supplementary Note 5). The computed REPs, which are normalized to the area of the measured traces atT = 12 T vib , are depicted as dashed lines in Fig. 5 and explicitly take into account for the participation of the entire multimode subspace to the Raman response 9 , with the adiabatic approximation that is limited to only one dimension, i.e., the one of the considered vibration. Building on this approach, when the mixing between ground and excited state normal modes is negligible (as confirmed in Supplementary Note 4 by evaluating the Duschinsky mixing 61 ), the unique free parameter required to extract the REPs is the displacement component along the GS normal mode under consideration which we obtained, along with the sign of the molecular polarizability derivative in Eq. (3), by time-dependent (TD) density functional theory (DFT) calculations (details are reported in the Methods section).
As shown in Fig. 5e, f, the normalized profiles and, most importantly, the sign of the measured traces (continuous curves) agree with the modeled signals (dashed lines), validating hence the sign of the reconstructed normal mode ES displacements. Small discrepancies between their lineshapes indicate a possible role of (1) the probe spectral  Finally, in order to evaluate the computed displacements, in Fig. 3e we compared the experimental absorption spectrum with the one obtained from the TD-DFT displacements (blue area and red line, respectively). Notably, the theoretical curve does not accurately reproduce the absorption profile, which in turn can be reproduced by proper scaling the low-frequency mode (ω g 0 g < 1000 cm −1 ) displacements by a factor 1.2 and the high-frequency ones by 1.4 (blue dashed line in Fig. 3e). In this respect, we note that the relative magnitude of the displacements can be directly extracted from the intensity of the ISRS oscillations, in view of their linear dependence on |d j |. Furthermore, their absolute amplitudes are accessed by taking advantage of the broadband nature of the ISRS detection, which enables a careful mapping of the REPs. The spectral profiles of these latter strongly depend on d j (details are reported in Supplementary Note 5), giving direct access to the displacement magnitudes.
We further note that the signal-to-noise ratio of the reconstructed REPs can be significantly improved by augmenting the statistical reliability of the extracted signal. To this aim, we reconstruct the REPs from multiple time delays (ΔT 2 400, 700 ½ fs with a 4 fs time-step) and PP chirps (C 2 = −4, 7, 30, 44, 71 fs 2 ). This can be done by regarding Eq. (6) (Methods section), for each fixed spectral component ω of the PP, as parametric in ΔT and C 2 , with the four unknowns <½R S g 0 ðωÞ, <½R AS g 0 ðωÞ, =½R S g 0 ðωÞ, =½R AS g 0 ðωÞ. Such over-determined inhomogeneous linear system can be simply solved in a least-squares sense, obtaining the best estimates for the four complex REP terms.
Absolute excited state geometries By using the procedure described above and taking into account for the 1-4 corrections, the complex resonant Raman excitation profiles of the five normal modes under consideration have been extracted and are reported in Fig. 6a. The corresponding excited state displacements with respect to the GS geometry are sketched in 6b and quantitatively reported in Table 2. Since the anti-Stokes REPs are shifted by one vibrational quantum with respect to the Stokes Raman excitation profiles (R S g 0 ðωÞ = R AS g 0 ðω À ω g 0 g Þ, see Eq. (4)), the former ones have been horizontally shifted for a more direct comparison of the measured traces with the REPs modeled from the absorption spectrum by using the transform method 9,60 and the measured displacements. The experimental traces are in are in good agreement with the theoretical profiles of all the normal modes under investigation, confirming the capability of the presented approach to fully capture the complex REPs and extract the corresponding displacements. Interestingly, our results reveal an increase of the distance between the oxygen and the distal carbon in the RhB central pyran ring, indicating an elongation of this latter in the excited state. Furthermore, the motion of the oxygen nuclei drives an out-of-plane rotation of the carboxylic acid, while a rearrangement of the lateral carbon atoms in the benzene rings leads to an outstretch expansion of these latters. Crucially, as anticipated in the introduction, by fully accessing the real and imaginary part of complex REPs, the reconstructed excitation profiles uniquely determine the sign of excited state displacements along the considered GS eigenvectors Q j . This ultimately unveils the molecular rearrangements in the electronic excited state along the five investigated normal coordinates, which are shown in Fig. 6b.  Table 1), are reported (as continuous lines) in e, f and compared with the REPs modeled by TD-DFT (dashed lines), normalized to the area of the measured traces atT = 12 T vib . Notably, the experimental and the theoretical signals share the same sign, validating hence the sign of the reconstructed normal mode ES displacement, while discrepancies between their spectral profiles suggest a role of the sample absorption, third-order dispersion and probe spectrum, which have not been included in the data analysis reported in this figure.
In summary, building on a diagrammatic treatment of the signal, we have demonstrated that the ISRS response encodes detailed information on both the real and the imaginary parts of molecular Raman excitation profiles. Based on a two-color ISRS scheme with an off-resonant Raman pump and a resonant probe pulse, we have shown how to experimentally measure complex REPs directly in the timedomain, determining the nuclear rearrangements on excited state potential energy surfaces. Notably, since the complex REP of a given normal mode is uniquely determined by two probe chirps and few time delays, the presented approach ensures fast acquisition times and high signal to noise ratios.
The proposed experimental scheme and theoretical model have been benchmarked in order to determine the ES of Rhodamine B dissolved in methanol projected onto different normal modes, identifying an elongation of the central pyran ring in the excited state, driven by an increased the distance between the oxygen and the distal carbon. We anticipate that applying the present protocol to access the ES displacements of photo-active compounds holds the potential for mapping in detail their reaction coordinates, unveiling the first stages of their relaxation dynamics.

Experimental setup
The Raman and the probe pulses used for the broadband ISRS experiment are synthesized from the same source, a Ti:sapphire laser, which generates 800 nm transform limited 40 fs pulses, with 1 kHz of repetition rate and an energy of 3.0 mJ. The RP is obtained by an inhouse built non-collinear optical parametric amplifier 63 , which generates vertically polarized pulses, centered at 700 nm. In order to avoid resonant excitation contributions from the pump, a long-pass filter is used to remove the RP spectral components below 600 nm, while a q (column 2) and in atomic units (column 3)  pair of chirped mirrors 64 is employed to control the RP compression, ensuring a 18-fs time duration. The vertically polarized broadband PP is synthesized from supercontinuum generation 45 by focusing a small portion of the laser source on a 3 mm sapphire crystal, producing a broadband white light continuum. An additional pair of chirped mirror is exploited to compress the generated PP, whose chirp is fine tuned by introducing thin glass windows on the probe optical path. The RP and the PP are focused in a non-collinear geometry (≈5°) on the rhodamine sample, which is contained in a 500 μm-thick quartz cuvette. The PP spectrum is monitored on a charge-coupled device upon frequency dispersion by a spectrometer, with a chopper synchronized at 500 Hz that blocks the Raman pump, in order to record consecutive probe pulse spectra, in the presence and in absence of the RP. Further details on the experimental scheme are reported in refs. 38, 52, 65.

Chirp dependence of the nonlinear ISRS response
The different contributions concurring to the ISRS signal can be computed via the Feynman diagrams introduced in Fig. 2. Specifically, they take into account for two classes of processes, depending on whether the two interactions with the off-resonant Raman pump (RP) involve the ket (Fig. 1a) or the bra-side of the density matrix (Fig. 1b). Among each class, the two interactions with the probe pulse (PP) project the system onto the entire electronic excited state manifold (with a transition amplitude ruled, according to the Franck-Condon principle 24 , by the overlap between the ground and excited states vibrational eigenfunctions). Importantly, the off-resonant RP ensures the generation of vibrational coherences jg 0 j ihg j j only on the initially populated electronic level for the A diagrams and jg j ihg 0 j j for the B, modulating the optical properties of the system at the frequencies ω j g 0 g of the stimulated normal modes ( j). In order to compactly evaluate the measured ISRS response, Eq.
Eq. (5) highlights that, while the vibrational excitations are triggered by the off-resonant Raman pump and are sensitive to the ground state molecular polarizability (via the ∂α ∂Q j term in Eq. (3)), the coherences are traced in a pump-probe fashion by the resonant probe pulse, which accesses the vibronic couplings via the Stokes and anti-Stokes Raman excitation profiles. Since in the considered experiments the RP and PP have parallel polarizations, the tensorial nature of the α term can be omitted, being the sample response sensitive only to the isotropic component of the polarizability.
To evaluate the ISRS dependence on the PP chirp 52,55,56,66,67 , the probe spectral field can be expressed as E P ðωÞ = E 0 P ðωÞe iC 2 ωÀω P ð Þ 2 , where E 0 P ðωÞ indicates the square root of the PP spectrum I P (ω). The ISRS signal can be conveniently recast by introducing the phase term ϕ = ϕðω, ΔT, C 2 Þ = ω g 0 g ΔT + 2 C 2 ðω À ω P Þω g 0 g =T Á ω g 0 g which corresponds to the product between the arrival time of the monitored probe wavelength, i.e.,T = ΔT + 2 C 2 ðω À ω P Þ and the frequency of the considered vibrational mode ω g 0 g . For a constant probe field (impulsive limit) in the region of interest (E 0 P ðωÞ = E 0 P ), it reads as (the complete derivation is reported in the Supplementary Note 1) This expression clarifies that tuning the probe chirp C 2 to selected values, slices of the dechirped signal (i.e., the ISRS response evaluated as a function of the time delay upon dechirpT) can be directly exploited to measure the sum and the difference of the real/imaginary parts of the Stokes/anti-Stokes complex REPs. We stress that, at odd with the spontaneous Raman case, the relative amplitudes of the impulsive Stokes/anti-Stokes Raman responses in Eqs. (2) and (6) do not depend on Boltzmann population factors, since vibrational excitations are stimulated by the off-resonant Raman pump and not by thermal excitation 68 .

Density functional theory calculations
DFT and TD-DFT calculations have been performed with CAM-B3LYP 69 functional and 6-311++g(2d,2p) basis set, taking into account for solvent effect by the polarizable continuum model, by using the Gaussian 16 software package 70 . Upon separate optimization of the ground and the excited state geometries by DFT/TD-DFT, the normal modes Q j and the associated transition dipole moments have been computed. Separate scans of the molecular geometries along the different GS eigenvectors Q j are performed to evaluate each molecular polarizability derivative ( ∂α ∂Q j ) in the considered reference frame. In parallel, the vertical projection of the ground state PES is monitored exploiting the ground state calculations to generate displaced geometries along the single normal modes under consideration. Further details are provided in the Supplementary Note 8.

Data availability
All the relevant data are available from the corresponding authors.